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A cluster algorithm is presented for the simulation of 
the g-state Potts models in which the number of spins is 
conserved in each state. The algorithm constructs Fortuin- 
Kasteleyn cluster configurations from spin configurations, in 
a way identical to the Swendsen-Wang algorithm; the spin as- 
signment to these clusters is however different, and conserves 
the number of spins for each state. Compared to traditional 
non-local spin-exchange algorithms, the cluster algorithm pre- 
sented here suffers less from critical slowing down, and con- 
sequently is more efficient near the critical temperature. 



I. INTRODUCTION 

Before 1987, the Potts model was almost exclusively 
simulated by means of the Metropolis algorithmtl, in 
which single-spin updates are proposed and either ac- 
cepted or rejected depending on the change in energy. 
This algorithm works quite satisfactorily, except close to 
the critical point. At the critical temperature, the cor- 
relation times increase with system size as L^, with a 
critical dynamic exponent equal to or slightly above two: 
for the two-states Potts model (Ising model), the critical 
dynamical exponent is reported to be z = 2.167 ± 0.001 
in two, and z = 2.02 ± 0.02 in three dimensionsH. The 
introduction of cluster algorithms has greatly advanced 
the accuracy with which critical properties of the Potts 
model and many other models in statistical physics can 
be studied. The first widely used cluster algorithm was 
introduced by Swendsen a nd W angBo; we will describe 
their algorithm in section II A . The dynamic exponent 
of the Swendsen-Wang algorithm in the two- and three- 
dimensional Ising model is reported to be z = 0.25 ±0.01 
and z = 0.54 ±0.02, respectively:B cluster algorithms are 
able to significantly reduce critical slowing down. 

To study multi-component lattice gases in the co- 
existence regime, for instance to study interfaces or equi- 
librium crystal shapes, one has to fix the number of par- 
ticles in the lattice gas for each component, i.e., the spin 
density for each state. One typically resorts to spin- 
exchange dynamics, with the unfortunate consequence 
of a critical slowing down at least as severe as experi- 
enced with the Metropolis algorithm applied to the regu- 
lar Potts model. The usual cluster algorithms do not 
conserve the spin densities. For the coasjsrved-order- 
parameter Ising model, Heringa and BlotaS'Ll recently in- 
troduced a cluster algorithm, which is in spirit related 



to the Wolff algorithnJa'D. It is reported to have hardly 
any critical slowing down, with a dynamical exponent of 
z = 0.21. This algorithm has not been generalized to 
Potts models with more than two states. 

In this paper, we present a modification of the 
Swendsen-Wang algorithm, to conserve the spin densi- 
ties. In the coming section, we describe their algorithm, 
and introduce our modified density-conserving cluster al- 
gorithm. In the next section, we present measurements 
of the critical dynamic exponent for our algorithm, and 
show its efficiency. The paper is concluded with a sum- 
mary and conclusions, and a discussion of future work. 



II. CLUSTER ALGORITHMS 

A. The Swendsen-Wang algorithm 

The Swendsen-Wang algorithm is designed to simulate 
the Potts model, defined by the Hamiltonian 
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in which J is the coupling constant, 5 denotes the Kro- 
necker delta function, and the summation runs over all 
pairs of nearest-neighbor sites, each having a spin with 
value [a = 1...Q). We use the usual symbols N for 
the number of lattice sites, L for the lateral dimension 
of the lattice with periodic boundary conditions, and 
Pi = (1/-/V) {ai, k) for the density of spins with value 
i. 

In this algorithm, the entire lattice is divided into clus- 
ters of aligned spins, to each of which a random new value 
is assigned. In detail, one step of the algorithm proceeds 
as follows: 

1. Visit all nearest-neighbor pairs of lattice sites; do 
nothing if the two spins are not aligned, but if they 
are, activate the bond between those two sites with 
a probabihty = 1 — exp(— /3J), where (3 is the 
inverse temperature. 

2. Group lattice sites that are connected by such ac- 
tivated bonds into clusters. 

3. Select a random new spin value for each cluster, 
and assign this spin value to each of the sites con- 
stituting the cluster. 

Steps 1, 2 and 3 are to be repeated many times, to obtain 
a set of sample configurations. 
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The proof of correctness for our density-conserving 
cluster algorithm is based on that for the Swendsen-Wang 
algorithm, which is presented in the remainder of this 
section. First we show detailed balance, next we discuss 
ergodicity. 

Suppose we denote the spin configuration before and 
after the move by Ca and Cb, respectively, with to- 
tal energies Ea and Ei,, and the intermediate clus- 
ter configuration Cm (also known as Fortuin-Kasteleyn 
representationa) ; furthermore, we write the probability 
to move from a configuration X to configuration Y as 
T{X Y). Then, the probability to move from a spin 
configuration Ca to a cluster configuration Cm is a prod- 
uct with factors Pc over all nearest-neighbor pairs of spins 
that are connected, times a product with factors 1 — Pc 
over all aligned nearest-neighbor pairs of spins that are 
disconnected: 
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and a similar expression for T{Cb Cm)- Since spins 
that are connected, are necessarily aligned both before 
and after the move, the first product on the right hand 
side is equal in T{Ca Cm) and T{Cb Cm)- AH 
factors in the second product on the right hand side 
dealing with pairs of spins that are aligned both before 
and after the move are also equal in T{Ca — > Cm) and 
T{Cb — > Cm)- That leaves in the ratio of the transition 
rates only the factors dealing with disconnected pairs of 
spins that are aligned either in configuration Ca, or in 
configuration Cb, but not both. The ratio of the tran- 
sition rates T{Ca Cm) and T{Cb — > Cm) therefore 
reduces to 

T{Ca > Cm) _ 



T{Cb C„ 
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Using that log(l — Pc) = —PJ, in combination with some 
rewriting, we obtain for the logarithm of this ratio 



l0g(T(Ca ^ Cm)) - l0g(r(Cb ^ Cm)) 
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As can easily been seen from the Hamiltonian eq. (|^), 
this is equal to — /3(i?a — Eb)- Since T{Cm Ca) = 
T{Cm Cb) — 2~" where n is the number of clusters in 
Cm, detailed balance follows: 



T{Cb — > Ca) _ T{Cb Cm) ■ T{Cm ^ Ca) 



T{Ca ^ Cb) T{Ca ^ Cm) ' T{C„ 

= exp{-(3{Ea-Eb))- 
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In addition to obeying detailed balance, the algorithm is 
ergodic, since there is a finite probability that in a given 
move all clusters will contain one site only, to which any 
value can be assigned. Since this algorithm is ergodic 
and satisfies detailed balance, it is guaranteed that even- 
tually these sample configurations are drawn from the 
Boltzmann distribution for the regular Potts model. The 
densities pi are not conserved in the Swendsen-Wang al- 
gorithm. 



B. Density-conserving cluster algorithm 

The topic of this paper is to present a modification 
to this algorithm, that ensures the conservation of the 
densities. This modification is made in step 3, in which 
the new spin values are assigned: rather than assigning 
random spin values to each cluster, we redistribute spin 
values over the clusters while conserving the spin densi- 
ties. As for the original Swendsen-Wang algorithm, the 
general idea is a two-step approach, Ca — > Cm — * Cb, 
where all the energetics required for obtaining detailed 
balance are incorporated in the construction of the clus- 
ters, and detailed balance is achieved by conservation of 
the property T(C„ ^ Ca) = T{Cm ^ Cb)- 

The first step towards such an algorithm is to devise 
an elementary move. The move we are looking for, is 
identifying one set of aligned clusters with spin value qi 
and another such set with spin value q2 ^ qi with exactly 
the same area (number of sites), and then exchanging the 
spin values qi and q2. 

How do we identify such sets? First of all, for each 
spin value i — 1 ... Q we group all clusters with spin 
value i into the set Si. Next, within each set we list 
these clusters in a random order, and keep track of the 
cumulative area. Every time that in two sets the same 
value for the cumulative area occurs, we have found an 
exchange point. If the spin values are exchanged in all 
clusters up to the exchange point, while the original spin 
values in all other clusters are conserved, the spin values 
of two sets of clusters are exchanged without violation of 
the spin density conservation. 

Unless extra measurements are taken, an algorithm 
based on these elementary moves will not obey detailed 
balance: the probability of occurrence for an exchange 
point is not necessarily equal before and after the cluster 
exchange. We denote the total number of clusters with 
spin- value q before the exchange takes place as Uq. Sup- 
pose that the exchange takes place between clusters with 
spin 1 and 2, and that the number of clusters with spin 
1, 2 that are to be exchanged is ai and 02, respectively, 
while the number of clusters with spin 1, 2 that are not 
to be exchanged is ni — ai and 712 — a2, respectively. 
The likelihood that there is an exchange point exactly 
between these sets of clusters is then equal to 
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while after the exchange, this probabihty becomes 
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To restore detailed balance, 
Metropolis acceptance ratio: 



it suffices to introduce a 



Pa = niin 
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Once this acceptance probability is included, the elemen- 
tary move can be used for a correct algorithm, since for 
two configurations X and Y , we now restored the prop- 
eviy T{Cm ^ X)^T{Cm-^Y). 

In an actual implementation, the total procedure is to 
make for each spin value a cumulative list of clusters, 
where the clusters are placed in a random order. Next, 
all exchange points are identified, the corresponding ex- 
changes are accepted with the probability as given in 
eq. (pf). It can be verified that also for the concatena- 
tion of exchange points, the product over all exchange 
points of the ratio of forward and backward acceptance 
probabilities, as given in eq. (||), equals 
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which exactly cancels the ratio of the number of ways 
in which the clusters can be sorted, i.e. the ratio of se- 
lection probabilities in forward and backward direction. 
Consequently T(C„. ^ Co) = T{Cm ^ Cb). 

The density-conserving algorithm is ergodic for the 
same reason that the Swendsen-Wang algorithm is er- 
godic: there is a finite probability that all clusters con- 
tain one site only, and then each of these can obtain any 
spin value (under the constraint on the densities). 

Having shown that the basic steps of our algorithm 
are correct we will now summarize the procedure in the 
form of a step-wise algorithm. Steps number 1 and 2 of 
the Swendsen-Wang algorithm remain unchanged. Step 
3 becomes: 

3a. For each state q, list all clusters with this spin- value 
in list 5*5, in a random order. 

3b. Order the lists with respect to the total area of their 
not-yet-assigned clusters. Use a random order for 
lists with equal such areas. If the first two lists 
(those with the largest and next-largest areas) are 
equal, exchange their colors with the probability as 
given by eq. (||) . Select one cluster from the first list 
and assign to it a new color. Update the ordering 
and repeat this step until spin values are assigned 
to all clusters. 



Note that the computational effort required for step 
3 scales with the total number of clusters n — '^^nq, 
whereas step 2 scales with the number of spins N in 
the system. Since n N, step 3 is repeated N/ (2n) 
times for each time step 2 is performed, and we still have 
an implementation in which the computational effort per 
sweep scales linearly with the number of sites; this greatly 
decreases the auto-correlation time. It actually also re- 
duces the dynamical critical exponent z, since the ratio 
N/n varies with the system size. 
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FIG. 1. Assignment of new spin values to the clusters 
in the three-states Potts model. The upper and lower part 
are the situations before and after cluster assignment, respec- 
tively. The length of each bar corresponds to the total mass 
PiN in the set 5*;. The different shades indicate different 
spin values before the assignment. The thick lines separate 
subsequent clusters in each list. The dotted lines indicate the 
exchange points A, B and C, where two of the masses coincide 
and after which the corresponding spin values are exchanged 
with the probability given in eq. (p|). 



III. COMPUTATIONAL PROPERTIES 

In order to compare the efficiency of the density- 
conserving cluster algorithm presented above with that of 
non-local spin-exchange (KawasakitH) dynamics, we have 
computed the energy autocorrelation times in the Ising 
model at critical temperature and equal spin densities, 
for several system sizes. Figures || and || show the corre- 
lation times as a function of the linear system size L of 
the two- and three-dimensional Ising model, respectively, 
both at their critical point. For all data points the corre- 
lation time r was obtained from a least-squares fit of the 
form e~*/'^, to the energy autocorrelation function. For 
the spin-exchange algorithm, these fits were done in the 
region where the autocorrelation drops from to e~^; 
for the cluster algorithm, we fitted in a broader region 
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(from 1 to e~"^), in order to have enough points to fit. All 
runs where started from a random configuration, which 
was thermalized over a time varying from 6000 MCS to 
60000 MCS. In order to generate enough statistics, the 
total length of the runs was set to ten times the thermal- 
ization time. The statistical errors were determined by 
repeating each run 10 to 50 times. 
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FIG. 2. correlation time r as a function of linear system 
size L for the two-dimensional Ising model, for spin-exchange 
dynamics (squares) and the magnetization-conserving cluster 
algorithm (circles). The lines have exponents of z = 2 and 
z = 0.38. 
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FIG. 3. correlation time r as a function of linear 
system size L for the three-dimensional Ising model, 
for spin-exchange dynamics (squares) and the magnetizar 
tion-conserving cluster algorithm (circles). The lines have 
exponents of 2; = 2 and z = 0.66. 



As expected, we find that the cluster algorithm suffers 
significantly less from critical slowing down and clearly 
outperforms spin-exchange dynamics at physically inter- 
esting lattice sizes in both two and three dimensions. 
Since one move in our cluster algorithm takes an amount 
of CPU-time comparable to what is required for one 
sweep in the non-local spin exchange, our cluster algo- 
rithm outperforms non-local spin exchange by one or two 
orders of magnitude, depending on the system size. 

For the non-local spin-exchange algorithm, we find 
a critical dynamic exponent of 2; = 2.0 in both two 
and three dimensions. This is in good agreement with 
the exponents of the three-dimensional non-conserving 
Metropolis algorithm (z 2.02 ± 0.02) but not with the 
critical exponent for the two-dimensional non- conserving 
Metropolis algorithm {z = 2.167 ± 0.001). Perhaps this 
is an indication that the conservation of the order pa- 
rameter affects the critical dynamical exponent, but our 
statistics are not conclusive. 

For the critical dynamic exponent for our new density- 
conserving cluster algorithm, we find values oi z = 0.38 ± 
0.01 in two, and z = 0.66 ± 0.02 in three dimensions. 
These values arc both slightly larger than non-conserved 
Swendsen-Wang values {z = 0.25 ± 0.01 and z = 0.54 ± 
0.02 for two and three dimensions respectively). 

IV. SUMMARY AND FUTURE WORK 

We have presented a density-conserving cluster algo- 
rithm for the Potts model. This algorithm is only mod- 
erately sensitive to critical slowing down: its dynamic 
critical exponent is found to be 2; = 0.38 ± 0.01 and 
z = 0.66 ± 0.02 for the two- and three-dimensional two- 
states Potts model, respectively. It outperforms the tra- 
ditional algorithm, non-local spin exchange, by one or 
two orders of magnitude. 

In future research, we will use this algorithm to study 
wetting properties, where the wetting takes place at a 
curved interface between two co-existing phases; such 
non-flat interfaces arise for instance between a droplet 
and a surrounding fluid. Other future applications will 
include the study of line tension between three co-existing 
phases, and equilibrium shapes in multi-component mix- 
tures. 
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